Calculate response (feeding ration, total weight of prey groups) and predictor variables for diet data, aggregate to get 1 stomach = 1 row
# Load libraries, install if needed
library(tidyverse); theme_set(theme_light(base_size = 10))
library(readxl)
library(tidylog)
library(RCurl)
library(viridis)
library(RColorBrewer)
library(patchwork)
library(janitor)
library(icesDatras)
library(mapdata)
library(patchwork)
library(rgdal)
library(raster)
library(sf)
library(rgeos)
library(chron)
library(lattice)
library(ncdf4)
library(marmap)
library(rnaturalearth)
library(rnaturalearthdata)
library(mapplots)
library(EnvStats)
library(qwraps2)
library(forcats)
#remotes::install_github("pbs-assess/sdmTMB")
library(sdmTMB)# To load entire cache in interactive r session, do: qwraps2::lazyload_cache_dir(path = "R/analysis/spatial_trend_models_cache/html")
# Specify map ranges
ymin = 55; ymax = 58; xmin = 14; xmax = 20
map_data <- rnaturalearth::ne_countries(
scale = "medium",
returnclass = "sf", continent = "europe")
# Crop the polygon for plotting and efficiency:
# st_bbox(map_data) # find the rough coordinates
swe_coast <- suppressWarnings(suppressMessages(
st_crop(map_data,
c(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax))))
# Transform our map into UTM 33 coordinates, which is the equal-area projection we fit in:
utm_zone33 <- 32633
swe_coast_proj <- sf::st_transform(swe_coast, crs = utm_zone33)
ggplot(swe_coast_proj) + geom_sf()
# Define plotting theme for main plot
theme_plot <- function(base_size = 10, base_family = "") {
theme_light(base_size = 10, base_family = "") +
theme(
axis.text.x = element_text(angle = 90),
axis.text = element_text(size = 8),
legend.text = element_text(size = 8),
legend.title = element_text(size = 8),
legend.position = "bottom",
legend.key.height = unit(0.2, "cm"),
legend.margin = margin(0, 0, 0, 0),
legend.box.margin = margin(-5, -5, -5, -5),
strip.text = element_text(size = 8, colour = 'black', margin = margin()),
strip.background = element_rect(fill = "grey90")
)
}
# Define plotting theme for facet_wrap map with years
theme_facet_map <- function(base_size = 10, base_family = "") {
theme_light(base_size = 10, base_family = "") +
theme(
axis.text.x = element_text(angle = 90),
axis.text = element_text(size = 6),
strip.text = element_text(size = 8, colour = 'black', margin = margin()),
strip.background = element_rect(fill = "grey90")
)
}
# Make default base map plot
plot_map_raster <-
ggplot(swe_coast_proj) +
geom_sf(size = 0.3) +
labs(x = "Longitude", y = "Latitude") +
theme_facet_map(base_size = 14)
# Function to convert from lat lon to UTM
LongLatToUTM <- function(x, y, zone){
xy <- data.frame(ID = 1:length(x), X = x, Y = y)
coordinates(xy) <- c("X", "Y")
proj4string(xy) <- CRS("+proj=longlat +datum=WGS84") ## for example
res <- spTransform(xy, CRS(paste("+proj=utm +zone=",zone," ellps=WGS84",sep='')))
return(as.data.frame(res))
}
d <- readr::read_csv("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/bentfish/data/for_analysis/stomach/full_stomach_data_21.10.25.csv") %>% dplyr::select(-X1)
#>
#> ── Column specification ────────────────────────────────────────────────────────
#> cols(
#> .default = col_double(),
#> Sample = col_character(),
#> Length.code = col_character(),
#> Prey.sp.code = col_character(),
#> Comment = col_character(),
#> Parasites.in.stomach = col_character(),
#> Coun. = col_character(),
#> Ship = col_logical(),
#> Cruise = col_character(),
#> Predator.code = col_character(),
#> Sex = col_character(),
#> Q.year = col_character(),
#> Date = col_date(format = ""),
#> Index = col_character(),
#> Ices.rect = col_character(),
#> source = col_character(),
#> Number = col_logical(),
#> Sample.type = col_logical(),
#> transect = col_logical(),
#> Depthstep = col_logical(),
#> Gonad.weight.roundfish = col_logical()
#> # ... with 7 more columns
#> )
#> ℹ Use `spec()` for the full column specifications.
# Drop NA coordinates (because these do not even have ices-rectangle, see diet data processing script in benthfish project)
d <- d %>% drop_na(Long) %>% drop_na(Lat)
#> drop_na: removed 853 rows (3%), 25,635 rows remaining
#> drop_na: no rows removed
# Add UTM coordinates (first rename )
utm_coords <- LongLatToUTM(d$Long, d$Lat, zone = 33)
d$X <- utm_coords$X/1000
d$Y <- utm_coords$Y/1000
head(data.frame(d))
#> SubDiv Year Month Day N.with.food N.regurgitated N.skeletal N.empty HN Sample
#> 1 28 2017 11 27 1 NA NA NA 52 507
#> 2 28 2017 11 27 1 NA NA NA 52 507
#> 3 28 2017 11 27 1 NA NA NA 52 507
#> 4 28 2017 11 27 1 NA NA NA 52 507
#> 5 28 2017 11 27 1 NA NA NA 52 507
#> 6 28 2017 11 27 1 NA NA NA 52 507
#> Length.code Prey.sp.code Prey.size Prey.weight Prey.nr Stage.digestion
#> 1 <NA> Saduria entomon 20 0.21 1 1
#> 2 <NA> Saduria entomon 16 0.10 1 1
#> 3 <NA> Saduria entomon 20 0.15 1 1
#> 4 <NA> Saduria entomon 20 0.17 1 1
#> 5 <NA> Saduria entomon NA 0.27 2 1
#> 6 <NA> Macoma balthica NA 0.19 1 1
#> Comment Parasites.in.stomach Quarter Coun. Ship Cruise Pred.size.mm
#> 1 <NA> <NA> 4 SWE NA BITS 200
#> 2 <NA> <NA> 4 SWE NA BITS 200
#> 3 <NA> <NA> 4 SWE NA BITS 200
#> 4 <NA> <NA> 4 SWE NA BITS 200
#> 5 <NA> <NA> 4 SWE NA BITS 200
#> 6 <NA> <NA> 4 SWE NA BITS 200
#> Pred.weight.g Predator.gutted.weight Predator.code Sex Maturity Gall.content
#> 1 100 92 FLE M 4 NA
#> 2 100 92 FLE M 4 NA
#> 3 100 92 FLE M 4 NA
#> 4 100 92 FLE M 4 NA
#> 5 100 92 FLE M 4 NA
#> 6 100 92 FLE M 4 NA
#> Q.year Age Date Index Lat Long Depth.catch Ices.rect
#> 1 4_2017 3 2017-11-27 OXBH.2017.52 57.46667 19.23333 67 43G9
#> 2 4_2017 3 2017-11-27 OXBH.2017.52 57.46667 19.23333 67 43G9
#> 3 4_2017 3 2017-11-27 OXBH.2017.52 57.46667 19.23333 67 43G9
#> 4 4_2017 3 2017-11-27 OXBH.2017.52 57.46667 19.23333 67 43G9
#> 5 4_2017 3 2017-11-27 OXBH.2017.52 57.46667 19.23333 67 43G9
#> 6 4_2017 3 2017-11-27 OXBH.2017.52 57.46667 19.23333 67 43G9
#> source Number Sample.type transect Depthstep Gonad.weight.roundfish
#> 1 Max Postdoc NA NA NA NA NA
#> 2 Max Postdoc NA NA NA NA NA
#> 3 Max Postdoc NA NA NA NA NA
#> 4 Max Postdoc NA NA NA NA NA
#> 5 Max Postdoc NA NA NA NA NA
#> 6 Max Postdoc NA NA NA NA NA
#> Liver.weight.roundfish Stomach.content Perc.stomac.content comments Processed
#> 1 NA NA NA NA NA
#> 2 NA NA NA NA NA
#> 3 NA NA NA NA NA
#> 4 NA NA NA NA NA
#> 5 NA NA NA NA NA
#> 6 NA NA NA NA NA
#> Unique.pred.id coord_from_ices_rect_centre X Y
#> 1 2017_4_52_507_FLE N 753.841 6377.247
#> 2 2017_4_52_507_FLE N 753.841 6377.247
#> 3 2017_4_52_507_FLE N 753.841 6377.247
#> 4 2017_4_52_507_FLE N 753.841 6377.247
#> 5 2017_4_52_507_FLE N 753.841 6377.247
#> 6 2017_4_52_507_FLE N 753.841 6377.247
sort(colnames(d))
#> [1] "Age" "Comment"
#> [3] "comments" "coord_from_ices_rect_centre"
#> [5] "Coun." "Cruise"
#> [7] "Date" "Day"
#> [9] "Depth.catch" "Depthstep"
#> [11] "Gall.content" "Gonad.weight.roundfish"
#> [13] "HN" "Ices.rect"
#> [15] "Index" "Lat"
#> [17] "Length.code" "Liver.weight.roundfish"
#> [19] "Long" "Maturity"
#> [21] "Month" "N.empty"
#> [23] "N.regurgitated" "N.skeletal"
#> [25] "N.with.food" "Number"
#> [27] "Parasites.in.stomach" "Perc.stomac.content"
#> [29] "Pred.size.mm" "Pred.weight.g"
#> [31] "Predator.code" "Predator.gutted.weight"
#> [33] "Prey.nr" "Prey.size"
#> [35] "Prey.sp.code" "Prey.weight"
#> [37] "Processed" "Q.year"
#> [39] "Quarter" "Sample"
#> [41] "Sample.type" "Sex"
#> [43] "Ship" "source"
#> [45] "Stage.digestion" "Stomach.content"
#> [47] "SubDiv" "transect"
#> [49] "Unique.pred.id" "X"
#> [51] "Y" "Year"
# [1] "Age" "Comment" "comments" "Coun."
# [5] "Cruise" "Date" "Day" "Depth.catch"
# [9] "Depthstep" "Gall.content" "Gonad.weight.roundfish" "HN"
# [13] "Ices.rect" "Index" "Lat" "Length.code"
# [17] "Liver.weight.roundfish" "Long" "Maturity" "Month"
# [21] "N.empty" "N.regurgitated" "N.skeletal" "N.with.food"
# [25] "Number" "Parasites.in.stomach" "Perc.stomac.content" "Pred.size.mm"
# [29] "Pred.weight.g" "Predator.code" "Predator.gutted.weight" "Prey.nr"
# [33] "Prey.size" "Prey.sp.code" "Prey.weight" "Processed"
# [37] "Q.year" "Quarter" "Sample" "Sample.type"
# [41] "Sex" "Ship" "source" "Stage.digestion"
# [45] "Stomach.content" "SubDiv" "transect" "Unique.pred.id"
# [49] "X"
# Plot data in space, color by survey
plot_map_raster +
geom_point(data = d, aes(x = X * 1000, y = Y * 1000, color = Cruise), size = 0.5) +
scale_color_brewer(palette = "Set2") +
facet_wrap(~ Cruise, ncol = 3) +
theme_plot()
# Calculate total weight of prey by predator ID and prey species (i.e., across prey sizes)
# Create wide data frame so that I can sum easily across prey groups (columns)
# I.e., one row = one stomach from here!
d_wide <- d %>%
drop_na(Prey.weight) %>%
group_by(Unique.pred.id, Prey.sp.code) %>%
summarise(tot_biom_per_prey = sum(Prey.weight)) %>%
ungroup() %>%
pivot_wider(names_from = Prey.sp.code, values_from = tot_biom_per_prey) %>%
mutate(across(everything(), ~replace_na(.x, 0))) %>% # Replace NA with 0, because it means the prey does not exist
janitor::clean_names()
#> drop_na: removed 1,216 rows (5%), 24,419 rows remaining
#> group_by: 2 grouping variables (Unique.pred.id, Prey.sp.code)
#> summarise: now 20,671 rows and 3 columns, one group variable remaining (Unique.pred.id)
#> ungroup: no grouping variables
#> pivot_wider: reorganized (Prey.sp.code, tot_biom_per_prey) into (Saduria entomon, Bylgides sarsi, Pontoporeia femorata, Clupea harengus, Mysis mixta, …) [was 20671x3, now 11189x117]
#> mutate: changed 8,723 values (78%) of 'Saduria entomon' (8723 fewer NA)
#> changed 8,675 values (78%) of 'Bylgides sarsi' (8675 fewer NA)
#> changed 10,425 values (93%) of 'Pontoporeia femorata' (10425 fewer NA)
#> changed 10,271 values (92%) of 'Clupea harengus' (10271 fewer NA)
#> changed 9,088 values (81%) of 'Mysis mixta' (9088 fewer NA)
#> changed 10,702 values (96%) of 'Gammarus sp.' (10702 fewer NA)
#> changed 9,454 values (84%) of 'Sprattus sprattus' (9454 fewer NA)
#> changed 11,151 values (>99%) of 'Enchelyopus cimbrius' (11151 fewer NA)
#> changed 10,857 values (97%) of 'Monoporeia affinis' (10857 fewer NA)
#> changed 11,184 values (>99%) of 'Mysis relicta' (11184 fewer NA)
#> changed 11,175 values (>99%) of 'Pontoporeia sp.' (11175 fewer NA)
#> changed 11,028 values (99%) of 'Gadus morhua' (11028 fewer NA)
#> changed 11,186 values (>99%) of 'Pisces Eggs' (11186 fewer NA)
#> changed 11,165 values (>99%) of 'Hyperia galba' (11165 fewer NA)
#> changed 10,336 values (92%) of 'Diastylis rathkei' (10336 fewer NA)
#> changed 10,486 values (94%) of 'Halicryptus spinulosus' (10486 fewer NA)
#> changed 10,184 values (91%) of 'Pisces' (10184 fewer NA)
#> changed 11,187 values (>99%) of 'Pomatoschistus minutus' (11187 fewer NA)
#> changed 10,206 values (91%) of 'Limecola balthica' (10206 fewer NA)
#> changed 11,186 values (>99%) of 'Mollusca' (11186 fewer NA)
#> changed 11,173 values (>99%) of 'Platichthys flesus' (11173 fewer NA)
#> changed 10,860 values (97%) of 'Gasterosteus aculeatus' (10860 fewer NA)
#> changed 10,875 values (97%) of 'Mysidae' (10875 fewer NA)
#> changed 10,985 values (98%) of 'Neomysis integer' (10985 fewer NA)
#> changed 10,723 values (96%) of 'Cumacea' (10723 fewer NA)
#> changed 11,185 values (>99%) of 'Corophium volutator' (11185 fewer NA)
#> changed 11,059 values (99%) of 'Stone' (11059 fewer NA)
#> changed 10,803 values (97%) of 'Clupeidae' (10803 fewer NA)
#> changed 10,777 values (96%) of 'Gobiidae' (10777 fewer NA)
#> changed 10,379 values (93%) of 'NA' (10379 fewer NA)
#> changed 11,143 values (>99%) of 'Scoloplos armiger' (11143 fewer NA)
#> changed 11,174 values (>99%) of 'plastic' (11174 fewer NA)
#> changed 11,094 values (99%) of 'Priapulus caudatus' (11094 fewer NA)
#> changed 11,129 values (99%) of 'Bivalvia' (11129 fewer NA)
#> changed 10,966 values (98%) of 'Mytilus sp.' (10966 fewer NA)
#> changed 10,892 values (97%) of 'Crangon crangon' (10892 fewer NA)
#> changed 11,140 values (>99%) of 'Idotea balthica' (11140 fewer NA)
#> changed 11,187 values (>99%) of 'Trachurus trachurus' (11187 fewer NA)
#> changed 11,187 values (>99%) of 'Annelida' (11187 fewer NA)
#> changed 11,095 values (99%) of 'Algae' (11095 fewer NA)
#> changed 11,185 values (>99%) of 'Priapulidae' (11185 fewer NA)
#> changed 11,182 values (>99%) of 'Waste' (11182 fewer NA)
#> changed 11,155 values (>99%) of 'Praunus flexuosus' (11155 fewer NA)
#> changed 11,178 values (>99%) of 'Unidentified mass' (11178 fewer NA)
#> changed 11,186 values (>99%) of 'Merlangius merlangus' (11186 fewer NA)
#> changed 11,122 values (99%) of 'Scales' (11122 fewer NA)
#> changed 11,186 values (>99%) of 'Gadidae' (11186 fewer NA)
#> changed 11,133 values (99%) of 'Crustacea' (11133 fewer NA)
#> changed 11,085 values (99%) of 'Amphipoda' (11085 fewer NA)
#> changed 11,135 values (>99%) of 'Spine' (11135 fewer NA)
#> changed 11,187 values (>99%) of 'Pleuronectes platessa' (11187 fewer NA)
#> changed 11,188 values (>99%) of 'Anguilla anguilla' (11188 fewer NA)
#> changed 11,188 values (>99%) of 'Empty' (11188 fewer NA)
#> changed 11,186 values (>99%) of 'Carbon' (11186 fewer NA)
#> changed 11,157 values (>99%) of 'Copepoda' (11157 fewer NA)
#> changed 11,180 values (>99%) of 'Zoarces viviparus' (11180 fewer NA)
#> changed 11,185 values (>99%) of 'Spinachia spinachia' (11185 fewer NA)
#> changed 11,186 values (>99%) of 'Pungitius pungitius' (11186 fewer NA)
#> changed 11,188 values (>99%) of 'Terebellides stroemii' (11188 fewer NA)
#> changed 11,171 values (>99%) of 'Palaemon sp.' (11171 fewer NA)
#> changed 11,183 values (>99%) of 'Ammodytes tobianus' (11183 fewer NA)
#> changed 11,187 values (>99%) of 'Cottidae' (11187 fewer NA)
#> changed 11,181 values (>99%) of 'Hediste diversicolor' (11181 fewer NA)
#> changed 11,174 values (>99%) of 'Palaemon elegans' (11174 fewer NA)
#> changed 11,164 values (>99%) of 'Ammodytidae' (11164 fewer NA)
#> changed 11,182 values (>99%) of 'Caridea' (11182 fewer NA)
#> changed 11,188 values (>99%) of 'Amphibalanus improvisus' (11188 fewer NA)
#> changed 11,188 values (>99%) of 'Myoxocephalus quadricornis' (11188 fewer NA)
#> changed 11,183 values (>99%) of 'Phyllodocida' (11183 fewer NA)
#> changed 11,104 values (99%) of 'Remains' (11104 fewer NA)
#> changed 11,180 values (>99%) of 'Palaemonidae' (11180 fewer NA)
#> changed 11,176 values (>99%) of 'Carcinus maenas' (11176 fewer NA)
#> changed 11,188 values (>99%) of 'Gastropoda' (11188 fewer NA)
#> changed 11,178 values (>99%) of 'Sand' (11178 fewer NA)
#> changed 11,186 values (>99%) of 'Cerastoderma glaucum' (11186 fewer NA)
#> changed 11,180 values (>99%) of 'Hyperoplus lanceolatus' (11180 fewer NA)
#> changed 11,163 values (>99%) of 'Mya arenaria' (11163 fewer NA)
#> changed 11,176 values (>99%) of 'Hydrobia sp.' (11176 fewer NA)
#> changed 11,182 values (>99%) of 'Wood' (11182 fewer NA)
#> changed 11,186 values (>99%) of 'Pleuronectiformes' (11186 fewer NA)
#> changed 11,188 values (>99%) of 'Scophthalmus maximus' (11188 fewer NA)
#> changed 11,161 values (>99%) of 'Polychaeta' (11161 fewer NA)
#> changed 11,091 values (99%) of 'Priapulida' (11091 fewer NA)
#> changed 11,188 values (>99%) of 'Halicryptus' (11188 fewer NA)
#> changed 11,187 values (>99%) of 'Neogobius melanostomus' (11187 fewer NA)
#> changed 11,187 values (>99%) of 'Gobius niger' (11187 fewer NA)
#> changed 11,188 values (>99%) of 'digestive tract' (11188 fewer NA)
#> changed 11,188 values (>99%) of 'Pleuronectidae' (11188 fewer NA)
#> changed 11,188 values (>99%) of 'sprattus sprattus' (11188 fewer NA)
#> changed 11,188 values (>99%) of 'Gasterosteidae' (11188 fewer NA)
#> changed 11,186 values (>99%) of 'litter' (11186 fewer NA)
#> changed 11,167 values (>99%) of 'Mysida' (11167 fewer NA)
#> changed 11,187 values (>99%) of 'Calanoida' (11187 fewer NA)
#> changed 11,188 values (>99%) of 'Belone belone' (11188 fewer NA)
#> changed 11,090 values (99%) of 'Pontoporeiidae' (11090 fewer NA)
#> changed 11,185 values (>99%) of 'Mucus' (11185 fewer NA)
#> changed 11,188 values (>99%) of 'Pectinaria sp.' (11188 fewer NA)
#> changed 11,069 values (99%) of 'Macoma balthica' (11069 fewer NA)
#> changed 10,904 values (97%) of 'remains' (10904 fewer NA)
#> changed 11,136 values (>99%) of 'stone' (11136 fewer NA)
#> changed 11,188 values (>99%) of 'carbon' (11188 fewer NA)
#> changed 11,187 values (>99%) of 'Nephtys ciliata' (11187 fewer NA)
#> changed 11,188 values (>99%) of 'Gastrosacus' (11188 fewer NA)
#> changed 11,188 values (>99%) of 'wood' (11188 fewer NA)
#> changed 11,188 values (>99%) of 'Litter/plastics' (11188 fewer NA)
#> changed 11,188 values (>99%) of 'clupeidae' (11188 fewer NA)
#> changed 11,176 values (>99%) of 'Pontoporeidae' (11176 fewer NA)
#> changed 11,182 values (>99%) of 'sand' (11182 fewer NA)
#> changed 11,187 values (>99%) of 'Decapoda' (11187 fewer NA)
#> changed 11,188 values (>99%) of 'Agonus cataphractus' (11188 fewer NA)
#> changed 11,188 values (>99%) of 'clupea harengus' (11188 fewer NA)
#> changed 11,173 values (>99%) of 'Halicryptus spinolusus' (11173 fewer NA)
#> changed 11,034 values (99%) of 'Mytilus edulis' (11034 fewer NA)
#> changed 11,181 values (>99%) of 'priapulida' (11181 fewer NA)
#> changed 11,188 values (>99%) of 'Prapulida' (11188 fewer NA)
#> changed 11,188 values (>99%) of 'Myoxocephalus scorpius' (11188 fewer NA)
head(d_wide)
#> # A tibble: 6 x 117
#> unique_pred_id saduria_entomon bylgides_sarsi pontoporeia_fem… clupea_harengus
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 1963_1_3_1963… 18.3 0 0 0
#> 2 1963_1_3_1963… 5.62 1.2 0 0
#> 3 1963_1_3_1963… 3.09 0 0 0
#> 4 1963_1_3_1963… 12.6 0 0 0
#> 5 1963_1_3_1963… 2.32 0 0 0
#> 6 1963_1_3_1963… 1.5 0 0 0
#> # … with 112 more variables: mysis_mixta <dbl>, gammarus_sp <dbl>,
#> # sprattus_sprattus <dbl>, enchelyopus_cimbrius <dbl>,
#> # monoporeia_affinis <dbl>, mysis_relicta <dbl>, pontoporeia_sp <dbl>,
#> # gadus_morhua <dbl>, pisces_eggs <dbl>, hyperia_galba <dbl>,
#> # diastylis_rathkei <dbl>, halicryptus_spinulosus <dbl>, pisces <dbl>,
#> # pomatoschistus_minutus <dbl>, limecola_balthica <dbl>, mollusca <dbl>,
#> # platichthys_flesus <dbl>, gasterosteus_aculeatus <dbl>, mysidae <dbl>,
#> # neomysis_integer <dbl>, cumacea <dbl>, corophium_volutator <dbl>,
#> # stone <dbl>, clupeidae <dbl>, gobiidae <dbl>, na <dbl>,
#> # scoloplos_armiger <dbl>, plastic <dbl>, priapulus_caudatus <dbl>,
#> # bivalvia <dbl>, mytilus_sp <dbl>, crangon_crangon <dbl>,
#> # idotea_balthica <dbl>, trachurus_trachurus <dbl>, annelida <dbl>,
#> # algae <dbl>, priapulidae <dbl>, waste <dbl>, praunus_flexuosus <dbl>,
#> # unidentified_mass <dbl>, merlangius_merlangus <dbl>, scales <dbl>,
#> # gadidae <dbl>, crustacea <dbl>, amphipoda <dbl>, spine <dbl>,
#> # pleuronectes_platessa <dbl>, anguilla_anguilla <dbl>, empty <dbl>,
#> # carbon <dbl>, copepoda <dbl>, zoarces_viviparus <dbl>,
#> # spinachia_spinachia <dbl>, pungitius_pungitius <dbl>,
#> # terebellides_stroemii <dbl>, palaemon_sp <dbl>, ammodytes_tobianus <dbl>,
#> # cottidae <dbl>, hediste_diversicolor <dbl>, palaemon_elegans <dbl>,
#> # ammodytidae <dbl>, caridea <dbl>, amphibalanus_improvisus <dbl>,
#> # myoxocephalus_quadricornis <dbl>, phyllodocida <dbl>, remains <dbl>,
#> # palaemonidae <dbl>, carcinus_maenas <dbl>, gastropoda <dbl>, sand <dbl>,
#> # cerastoderma_glaucum <dbl>, hyperoplus_lanceolatus <dbl>,
#> # mya_arenaria <dbl>, hydrobia_sp <dbl>, wood <dbl>, pleuronectiformes <dbl>,
#> # scophthalmus_maximus <dbl>, polychaeta <dbl>, priapulida <dbl>,
#> # halicryptus <dbl>, neogobius_melanostomus <dbl>, gobius_niger <dbl>,
#> # digestive_tract <dbl>, pleuronectidae <dbl>, sprattus_sprattus_2 <dbl>,
#> # gasterosteidae <dbl>, litter <dbl>, mysida <dbl>, calanoida <dbl>,
#> # belone_belone <dbl>, pontoporeiidae <dbl>, mucus <dbl>,
#> # pectinaria_sp <dbl>, macoma_balthica <dbl>, remains_2 <dbl>, stone_2 <dbl>,
#> # carbon_2 <dbl>, nephtys_ciliata <dbl>, gastrosacus <dbl>, wood_2 <dbl>, …
str(d_wide)
#> tibble [11,189 × 117] (S3: tbl_df/tbl/data.frame)
#> $ unique_pred_id : chr [1:11189] "1963_1_3_1963_3_18_3_MAZ_41G8_22_NA_500_NA_1_COD" "1963_1_3_1963_3_18_3_MAZ_41G8_39_NA_450_NA_1_COD" "1963_1_3_1963_3_18_3_MAZ_41G8_41_NA_560_NA_1_COD" "1963_1_3_1963_3_18_3_MAZ_41G8_42_NA_390_NA_1_COD" ...
#> $ saduria_entomon : num [1:11189] 18.27 5.62 3.09 12.6 2.32 ...
#> $ bylgides_sarsi : num [1:11189] 0 1.2 0 0 0 0 0 0 0 0 ...
#> $ pontoporeia_femorata : num [1:11189] 0 0 0 0 0 0 3.34 0 0 0.1 ...
#> $ clupea_harengus : num [1:11189] 0 0 0 0 0 ...
#> $ mysis_mixta : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ gammarus_sp : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ sprattus_sprattus : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ enchelyopus_cimbrius : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ monoporeia_affinis : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ mysis_relicta : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ pontoporeia_sp : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ gadus_morhua : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ pisces_eggs : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ hyperia_galba : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ diastylis_rathkei : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ halicryptus_spinulosus : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ pisces : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ pomatoschistus_minutus : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ limecola_balthica : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ mollusca : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ platichthys_flesus : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ gasterosteus_aculeatus : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ mysidae : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ neomysis_integer : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ cumacea : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ corophium_volutator : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ stone : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ clupeidae : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ gobiidae : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ na : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ scoloplos_armiger : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ plastic : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ priapulus_caudatus : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ bivalvia : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ mytilus_sp : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ crangon_crangon : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ idotea_balthica : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ trachurus_trachurus : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ annelida : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ algae : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ priapulidae : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ waste : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ praunus_flexuosus : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ unidentified_mass : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ merlangius_merlangus : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ scales : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ gadidae : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ crustacea : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ amphipoda : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ spine : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ pleuronectes_platessa : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ anguilla_anguilla : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ empty : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ carbon : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ copepoda : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ zoarces_viviparus : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ spinachia_spinachia : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ pungitius_pungitius : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ terebellides_stroemii : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ palaemon_sp : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ ammodytes_tobianus : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ cottidae : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ hediste_diversicolor : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ palaemon_elegans : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ ammodytidae : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ caridea : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ amphibalanus_improvisus : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ myoxocephalus_quadricornis: num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ phyllodocida : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ remains : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ palaemonidae : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ carcinus_maenas : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ gastropoda : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ sand : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ cerastoderma_glaucum : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ hyperoplus_lanceolatus : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ mya_arenaria : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ hydrobia_sp : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ wood : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ pleuronectiformes : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ scophthalmus_maximus : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ polychaeta : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ priapulida : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ halicryptus : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ neogobius_melanostomus : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ gobius_niger : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ digestive_tract : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ pleuronectidae : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ sprattus_sprattus_2 : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ gasterosteidae : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ litter : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ mysida : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ calanoida : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ belone_belone : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ pontoporeiidae : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ mucus : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ pectinaria_sp : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> $ macoma_balthica : num [1:11189] 0 0 0 0 0 0 0 0 0 0 ...
#> [list output truncated]
colnames(d_wide)
#> [1] "unique_pred_id" "saduria_entomon"
#> [3] "bylgides_sarsi" "pontoporeia_femorata"
#> [5] "clupea_harengus" "mysis_mixta"
#> [7] "gammarus_sp" "sprattus_sprattus"
#> [9] "enchelyopus_cimbrius" "monoporeia_affinis"
#> [11] "mysis_relicta" "pontoporeia_sp"
#> [13] "gadus_morhua" "pisces_eggs"
#> [15] "hyperia_galba" "diastylis_rathkei"
#> [17] "halicryptus_spinulosus" "pisces"
#> [19] "pomatoschistus_minutus" "limecola_balthica"
#> [21] "mollusca" "platichthys_flesus"
#> [23] "gasterosteus_aculeatus" "mysidae"
#> [25] "neomysis_integer" "cumacea"
#> [27] "corophium_volutator" "stone"
#> [29] "clupeidae" "gobiidae"
#> [31] "na" "scoloplos_armiger"
#> [33] "plastic" "priapulus_caudatus"
#> [35] "bivalvia" "mytilus_sp"
#> [37] "crangon_crangon" "idotea_balthica"
#> [39] "trachurus_trachurus" "annelida"
#> [41] "algae" "priapulidae"
#> [43] "waste" "praunus_flexuosus"
#> [45] "unidentified_mass" "merlangius_merlangus"
#> [47] "scales" "gadidae"
#> [49] "crustacea" "amphipoda"
#> [51] "spine" "pleuronectes_platessa"
#> [53] "anguilla_anguilla" "empty"
#> [55] "carbon" "copepoda"
#> [57] "zoarces_viviparus" "spinachia_spinachia"
#> [59] "pungitius_pungitius" "terebellides_stroemii"
#> [61] "palaemon_sp" "ammodytes_tobianus"
#> [63] "cottidae" "hediste_diversicolor"
#> [65] "palaemon_elegans" "ammodytidae"
#> [67] "caridea" "amphibalanus_improvisus"
#> [69] "myoxocephalus_quadricornis" "phyllodocida"
#> [71] "remains" "palaemonidae"
#> [73] "carcinus_maenas" "gastropoda"
#> [75] "sand" "cerastoderma_glaucum"
#> [77] "hyperoplus_lanceolatus" "mya_arenaria"
#> [79] "hydrobia_sp" "wood"
#> [81] "pleuronectiformes" "scophthalmus_maximus"
#> [83] "polychaeta" "priapulida"
#> [85] "halicryptus" "neogobius_melanostomus"
#> [87] "gobius_niger" "digestive_tract"
#> [89] "pleuronectidae" "sprattus_sprattus_2"
#> [91] "gasterosteidae" "litter"
#> [93] "mysida" "calanoida"
#> [95] "belone_belone" "pontoporeiidae"
#> [97] "mucus" "pectinaria_sp"
#> [99] "macoma_balthica" "remains_2"
#> [101] "stone_2" "carbon_2"
#> [103] "nephtys_ciliata" "gastrosacus"
#> [105] "wood_2" "litter_plastics"
#> [107] "clupeidae_2" "pontoporeidae"
#> [109] "sand_2" "decapoda"
#> [111] "agonus_cataphractus" "clupea_harengus_2"
#> [113] "halicryptus_spinolusus" "mytilus_edulis"
#> [115] "priapulida_2" "prapulida"
#> [117] "myoxocephalus_scorpius"
# Now make some calculations and aggregate some taxonomic level. Since all columns are assigned to
# some higher level group (or the same group), the sum of these is the total stomach content.
# Note that I have one group for unidentified clupeids, but also sprat and herring. So if I want the
# total of some aggregated group, then I need to add all the sub-groups.
d_wide2 <- d_wide %>%
mutate(amphipoda_tot = bylgides_sarsi + hyperia_galba + gammarus_sp + monoporeia_affinis +
corophium_volutator + amphipoda,
bivalvia_tot = bivalvia + mytilus_sp + cerastoderma_glaucum + mya_arenaria + macoma_balthica +
mytilus_edulis,
clupeidae_tot = clupeidae + clupeidae_2,
sprattus_sprattus_tot = sprattus_sprattus + sprattus_sprattus_2,
clupea_harengus_tot = clupea_harengus + clupea_harengus_2,
gadus_morhua_tot = gadus_morhua,
gadiformes_tot = gadidae + merlangius_merlangus,
gobiidae_tot = gobiidae,
mysidae_tot = mysidae + neomysis_integer + mysis_mixta + mysida + gastrosacus,
non_bio_tot = stone + plastic + sand + wood + litter + carbon + stone_2 + carbon_2 + wood_2 +
litter_plastics + sand_2,
other_crustacea_tot = pontoporeia_femorata + crangon_crangon + idotea_balthica +
praunus_flexuosus + crustacea + diastylis_rathkei + palaemon_sp + palaemon_elegans + caridea +
amphibalanus_improvisus + palaemonidae + carcinus_maenas + copepoda + calanoida + pontoporeiidae + decapoda,
other_tot = halicryptus_spinulosus + priapulus_caudatus + annelida + algae + priapulidae + waste +
unidentified_mass + spine + empty + mollusca + na + remains + gastropoda + hydrobia_sp +
priapulida + halicryptus + digestive_tract + mucus + remains_2 + pontoporeidae +
halicryptus_spinolusus + priapulida_2 + prapulida,
other_pisces_tot = enchelyopus_cimbrius + trachurus_trachurus + gasterosteus_aculeatus +
scales + pleuronectes_platessa + anguilla_anguilla + pungitius_pungitius + ammodytes_tobianus +
cottidae + spinachia_spinachia + zoarces_viviparus + ammodytidae + myoxocephalus_quadricornis +
hyperoplus_lanceolatus + pleuronectiformes + scophthalmus_maximus + neogobius_melanostomus + gobius_niger +
pleuronectidae + gasterosteidae + belone_belone + agonus_cataphractus + myoxocephalus_scorpius,
platichthys_flesus_tot = platichthys_flesus,
polychaeta_tot = bylgides_sarsi + scoloplos_armiger + terebellides_stroemii + hediste_diversicolor +
phyllodocida + polychaeta + pectinaria_sp + nephtys_ciliata,
saduria_entomon_tot = saduria_entomon
)
#> mutate: new variable 'amphipoda_tot' (double) with 784 unique values and 0% NA
#> new variable 'bivalvia_tot' (double) with 257 unique values and 0% NA
#> new variable 'clupeidae_tot' (double) with 297 unique values and 0% NA
#> new variable 'sprattus_sprattus_tot' (double) with 1,295 unique values and 0% NA
#> new variable 'clupea_harengus_tot' (double) with 826 unique values and 0% NA
#> new variable 'gadus_morhua_tot' (double) with 156 unique values and 0% NA
#> new variable 'gadiformes_tot' (double) with 7 unique values and 0% NA
#> new variable 'gobiidae_tot' (double) with 165 unique values and 0% NA
#> new variable 'mysidae_tot' (double) with 445 unique values and 0% NA
#> new variable 'non_bio_tot' (double) with 66 unique values and 0% NA
#> new variable 'other_crustacea_tot' (double) with 304 unique values and 0% NA
#> new variable 'other_tot' (double) with 304 unique values and 0% NA
#> new variable 'other_pisces_tot' (double) with 312 unique values and 0% NA
#> new variable 'platichthys_flesus_tot' (double) with 17 unique values and 0% NA
#> new variable 'polychaeta_tot' (double) with 738 unique values and 0% NA
#> new variable 'saduria_entomon_tot' (double) with 784 unique values and 0% NA
# 16 prey groups in total
length(unique(d$Unique.pred.id))
#> [1] 12310
nrow(d_wide2) # The reason they differ in length (nrow) is because I remove NA prey weight
#> [1] 11189
length(unique(drop_na(d, Prey.weight)$Unique.pred.id))
#> drop_na: removed 1,216 rows (5%), 24,419 rows remaining
#> [1] 11189
# Select only columns aggregated columns (ending with _tot)
data.frame(d_wide2[1, c(1, 118:133)])
#> unique_pred_id amphipoda_tot bivalvia_tot
#> 1 1963_1_3_1963_3_18_3_MAZ_41G8_22_NA_500_NA_1_COD 0 0
#> clupeidae_tot sprattus_sprattus_tot clupea_harengus_tot gadus_morhua_tot
#> 1 0 0 0 0
#> gadiformes_tot gobiidae_tot mysidae_tot non_bio_tot other_crustacea_tot
#> 1 0 0 0 0 0
#> other_tot other_pisces_tot platichthys_flesus_tot polychaeta_tot
#> 1 0 0 0 0
#> saduria_entomon_tot
#> 1 18.27
d_wide3 <- d_wide2 %>% dplyr::select(c(1, 118:133))
sort(colnames(d_wide3))
#> [1] "amphipoda_tot" "bivalvia_tot" "clupea_harengus_tot"
#> [4] "clupeidae_tot" "gadiformes_tot" "gadus_morhua_tot"
#> [7] "gobiidae_tot" "mysidae_tot" "non_bio_tot"
#> [10] "other_crustacea_tot" "other_pisces_tot" "other_tot"
#> [13] "platichthys_flesus_tot" "polychaeta_tot" "saduria_entomon_tot"
#> [16] "sprattus_sprattus_tot" "unique_pred_id"
# Add back in other information
d_sub <- d %>%
dplyr::select(Year, Quarter, Cruise, HN, Sample, Predator.code, X, Y, Lat, Long, Ices.rect,
Pred.size.mm, Pred.weight.g, Unique.pred.id, source) %>%
rename("Predator.spec" = "Predator.code") %>%
distinct(Unique.pred.id, .keep_all = TRUE) %>%
janitor::clean_names()
#> rename: renamed one variable (Predator.spec)
#> distinct: removed 13,325 rows (52%), 12,310 rows remaining
d_wide4 <- left_join(d_wide3, d_sub) # Why missing 1,121 IDs? I lose them when making d_wide, I filter non-NA prey weights!
#> Joining, by = "unique_pred_id"
#> left_join: added 14 columns (year, quarter, cruise, hn, sample, …)
#> > rows only in x 0
#> > rows only in y ( 1,121)
#> > matched rows 11,189
#> > ========
#> > rows total 11,189
# missing_ids <- unique(d_sub$unique_pred_id)[!unique(d_sub$unique_pred_id) %in% unique(d_wide3$unique_pred_id)]
# d %>% filter(Unique.pred.id %in% missing_ids) %>% distinct(Prey.weight)
colnames(d_wide4)
#> [1] "unique_pred_id" "amphipoda_tot" "bivalvia_tot"
#> [4] "clupeidae_tot" "sprattus_sprattus_tot" "clupea_harengus_tot"
#> [7] "gadus_morhua_tot" "gadiformes_tot" "gobiidae_tot"
#> [10] "mysidae_tot" "non_bio_tot" "other_crustacea_tot"
#> [13] "other_tot" "other_pisces_tot" "platichthys_flesus_tot"
#> [16] "polychaeta_tot" "saduria_entomon_tot" "year"
#> [19] "quarter" "cruise" "hn"
#> [22] "sample" "predator_spec" "x"
#> [25] "y" "lat" "long"
#> [28] "ices_rect" "pred_size_mm" "pred_weight_g"
#> [31] "source"
data.frame(d_wide4[1, 2:17])
#> amphipoda_tot bivalvia_tot clupeidae_tot sprattus_sprattus_tot
#> 1 0 0 0 0
#> clupea_harengus_tot gadus_morhua_tot gadiformes_tot gobiidae_tot mysidae_tot
#> 1 0 0 0 0 0
#> non_bio_tot other_crustacea_tot other_tot other_pisces_tot
#> 1 0 0 0 0
#> platichthys_flesus_tot polychaeta_tot saduria_entomon_tot
#> 1 0 0 18.27
d_wide4 <- d_wide4 %>% mutate(tot_prey_biom = rowSums(.[2:17]))
#> mutate: new variable 'tot_prey_biom' (double) with 3,159 unique values and 0% NA
d_wide4 %>% group_by(unique_pred_id) %>% mutate(n = n()) %>% ungroup() %>% distinct(n)
#> group_by: one grouping variable (unique_pred_id)
#> mutate (grouped): new variable 'n' (integer) with one unique value and 0% NA
#> ungroup: no grouping variables
#> distinct: removed 11,188 rows (>99%), one row remaining
#> # A tibble: 1 x 1
#> n
#> <int>
#> 1 1
# Split data by species
fle <- d_wide4 %>%
filter(predator_spec == "FLE") %>%
mutate(pred_cm = pred_size_mm/10,
pred_cm_class = cut(pred_cm, breaks = c(0, 9, 20, 30, 200), right = TRUE))
#> filter: removed 8,932 rows (80%), 2,257 rows remaining
#> mutate: new variable 'pred_cm' (double) with 113 unique values and <1% NA
#> new variable 'pred_cm_class' (factor) with 5 unique values and <1% NA
head(data.frame(fle), 20)
#> unique_pred_id amphipoda_tot bivalvia_tot clupeidae_tot
#> 1 2015_2_1010_4_FLE 0 2.60 0
#> 2 2015_2_1010_402_FLE 0 0.00 0
#> 3 2015_2_108_3_FLE 0 0.00 0
#> 4 2015_2_202_403_FLE 0 5.46 0
#> 5 2015_2_203_1_FLE 0 1.74 0
#> 6 2015_2_204_2_FLE 0 0.04 0
#> 7 2015_2_501_10_FLE 0 0.99 0
#> 8 2015_2_501_60_FLE 0 0.00 0
#> 9 2015_2_501_8_FLE 0 0.00 0
#> 10 2015_2_501_9_FLE 0 0.00 0
#> 11 2015_2_504_5_FLE 0 0.00 0
#> 12 2015_2_504_58_FLE 0 0.00 0
#> 13 2015_2_504_6_FLE 0 0.00 0
#> 14 2015_2_504_7_FLE 0 0.00 0
#> 15 2015_2_506_11_FLE 0 0.00 0
#> 16 2015_2_506_12_FLE 0 0.00 0
#> 17 2015_2_506_13_FLE 0 0.00 0
#> 18 2015_2_506_14_FLE 0 0.00 0
#> 19 2015_2_506_15_FLE 0 0.00 0
#> 20 2015_2_506_16_FLE 0 0.00 0
#> sprattus_sprattus_tot clupea_harengus_tot gadus_morhua_tot gadiformes_tot
#> 1 0 0 0 0
#> 2 0 0 0 0
#> 3 0 0 0 0
#> 4 0 0 0 0
#> 5 0 0 0 0
#> 6 0 0 0 0
#> 7 0 0 0 0
#> 8 0 0 0 0
#> 9 0 0 0 0
#> 10 0 0 0 0
#> 11 0 0 0 0
#> 12 0 0 0 0
#> 13 0 0 0 0
#> 14 0 0 0 0
#> 15 0 0 0 0
#> 16 0 0 0 0
#> 17 0 0 0 0
#> 18 0 0 0 0
#> 19 0 0 0 0
#> 20 0 0 0 0
#> gobiidae_tot mysidae_tot non_bio_tot other_crustacea_tot other_tot
#> 1 0 0 0.19 0.00 0.00
#> 2 0 0 0.00 0.00 0.00
#> 3 0 0 0.00 0.00 0.00
#> 4 0 0 0.13 0.00 0.00
#> 5 0 0 0.00 0.00 0.00
#> 6 0 0 0.00 0.00 0.00
#> 7 0 0 0.00 0.00 0.00
#> 8 0 0 0.00 0.00 0.00
#> 9 0 0 0.00 0.01 0.03
#> 10 0 0 0.00 0.00 0.00
#> 11 0 0 0.00 0.02 0.00
#> 12 0 0 0.00 0.00 0.02
#> 13 0 0 0.00 0.00 0.00
#> 14 0 0 0.00 0.00 0.00
#> 15 0 0 0.00 0.00 0.00
#> 16 0 0 0.00 0.10 0.00
#> 17 0 0 0.00 0.00 0.00
#> 18 0 0 0.00 0.00 0.00
#> 19 0 0 0.00 0.00 0.00
#> 20 0 0 0.00 0.00 0.00
#> other_pisces_tot platichthys_flesus_tot polychaeta_tot saduria_entomon_tot
#> 1 2.93 0 0 0.00
#> 2 0.00 0 0 0.00
#> 3 0.00 0 0 5.61
#> 4 0.00 0 0 1.02
#> 5 0.00 0 0 3.10
#> 6 0.00 0 0 2.98
#> 7 0.00 0 0 0.00
#> 8 0.00 0 0 0.00
#> 9 0.00 0 0 0.00
#> 10 0.00 0 0 0.00
#> 11 0.00 0 0 0.08
#> 12 0.00 0 0 0.00
#> 13 0.00 0 0 0.00
#> 14 0.00 0 0 0.00
#> 15 0.00 0 0 0.00
#> 16 0.00 0 0 0.00
#> 17 0.00 0 0 0.00
#> 18 0.00 0 0 0.00
#> 19 0.00 0 0 0.00
#> 20 0.00 0 0 0.00
#> year quarter cruise hn sample predator_spec x y lat long
#> 1 2015 2 INSPIRE 1010 4 FLE 486.2897 6209.440 56.03 14.78
#> 2 2015 2 INSPIRE 1010 402 FLE 486.2897 6209.440 56.03 14.78
#> 3 2015 2 INSPIRE 108 3 FLE 483.1607 6206.112 56.00 14.73
#> 4 2015 2 INSPIRE 202 403 FLE 485.6369 6200.539 55.95 14.77
#> 5 2015 2 INSPIRE 203 1 FLE 486.9129 6209.438 56.03 14.79
#> 6 2015 2 INSPIRE 204 2 FLE 485.6517 6204.990 55.99 14.77
#> 7 2015 2 INSPIRE 501 10 FLE 478.0250 6177.198 55.74 14.65
#> 8 2015 2 INSPIRE 501 60 FLE 478.0250 6177.198 55.74 14.65
#> 9 2015 2 INSPIRE 501 8 FLE 478.0250 6177.198 55.74 14.65
#> 10 2015 2 INSPIRE 501 9 FLE 478.0250 6177.198 55.74 14.65
#> 11 2015 2 INSPIRE 504 5 FLE 477.9969 6171.634 55.69 14.65
#> 12 2015 2 INSPIRE 504 58 FLE 477.9969 6171.634 55.69 14.65
#> 13 2015 2 INSPIRE 504 6 FLE 477.9969 6171.634 55.69 14.65
#> 14 2015 2 INSPIRE 504 7 FLE 477.9969 6171.634 55.69 14.65
#> 15 2015 2 INSPIRE 506 11 FLE 479.2754 6176.079 55.73 14.67
#> 16 2015 2 INSPIRE 506 12 FLE 479.2754 6176.079 55.73 14.67
#> 17 2015 2 INSPIRE 506 13 FLE 479.2754 6176.079 55.73 14.67
#> 18 2015 2 INSPIRE 506 14 FLE 479.2754 6176.079 55.73 14.67
#> 19 2015 2 INSPIRE 506 15 FLE 479.2754 6176.079 55.73 14.67
#> 20 2015 2 INSPIRE 506 16 FLE 479.2754 6176.079 55.73 14.67
#> ices_rect pred_size_mm pred_weight_g source tot_prey_biom pred_cm
#> 1 41G4 270 236 Kevin's MSC 5.72 27.0
#> 2 41G4 304 476 Kevin's MSC 0.00 30.4
#> 3 41G4 350 334 Kevin's MSC 5.61 35.0
#> 4 40G4 302 264 Kevin's MSC 6.61 30.2
#> 5 41G4 299 246 Kevin's MSC 4.84 29.9
#> 6 40G4 346 346 Kevin's MSC 3.02 34.6
#> 7 40G4 314 305 Kevin's MSC 0.99 31.4
#> 8 40G4 330 548 Kevin's MSC 0.00 33.0
#> 9 40G4 202 89 Kevin's MSC 0.04 20.2
#> 10 40G4 298 251 Kevin's MSC 0.00 29.8
#> 11 40G4 352 325 Kevin's MSC 0.10 35.2
#> 12 40G4 313 319 Kevin's MSC 0.02 31.3
#> 13 40G4 328 367 Kevin's MSC 0.00 32.8
#> 14 40G4 250 172 Kevin's MSC 0.00 25.0
#> 15 40G4 299 196 Kevin's MSC 0.00 29.9
#> 16 40G4 280 184 Kevin's MSC 0.10 28.0
#> 17 40G4 349 389 Kevin's MSC 0.00 34.9
#> 18 40G4 335 368 Kevin's MSC 0.00 33.5
#> 19 40G4 340 371 Kevin's MSC 0.00 34.0
#> 20 40G4 337 380 Kevin's MSC 0.00 33.7
#> pred_cm_class
#> 1 (20,30]
#> 2 (30,200]
#> 3 (30,200]
#> 4 (30,200]
#> 5 (20,30]
#> 6 (30,200]
#> 7 (30,200]
#> 8 (30,200]
#> 9 (20,30]
#> 10 (20,30]
#> 11 (30,200]
#> 12 (30,200]
#> 13 (30,200]
#> 14 (20,30]
#> 15 (20,30]
#> 16 (20,30]
#> 17 (30,200]
#> 18 (30,200]
#> 19 (30,200]
#> 20 (30,200]
unique(fle$pred_cm_class)
#> [1] (20,30] (30,200] (9,20] (0,9] <NA>
#> Levels: (0,9] (9,20] (20,30] (30,200]
cod <- d_wide4 %>%
filter(predator_spec == "COD") %>%
mutate(pred_cm = pred_size_mm/10,
pred_cm_class = cut(pred_cm, breaks = c(0, 6, 20, 30, 40, 50, 200), right = TRUE))
#> filter: removed 2,257 rows (20%), 8,932 rows remaining
#> mutate: new variable 'pred_cm' (double) with 224 unique values and <1% NA
#> new variable 'pred_cm_class' (factor) with 7 unique values and 3% NA
head(data.frame(cod), 20)
#> unique_pred_id amphipoda_tot bivalvia_tot
#> 1 1963_1_3_1963_3_18_3_MAZ_41G8_22_NA_500_NA_1_COD 0.0 0
#> 2 1963_1_3_1963_3_18_3_MAZ_41G8_39_NA_450_NA_1_COD 1.2 0
#> 3 1963_1_3_1963_3_18_3_MAZ_41G8_41_NA_560_NA_1_COD 0.0 0
#> 4 1963_1_3_1963_3_18_3_MAZ_41G8_42_NA_390_NA_1_COD 0.0 0
#> 5 1963_1_3_1963_3_18_3_MAZ_41G8_45_NA_270_NA_1_COD 0.0 0
#> 6 1963_1_3_1963_3_18_3_MAZ_41G8_47_NA_270_NA_1_COD 0.0 0
#> 7 1963_1_3_1963_3_18_3_MAZ_41G8_51_NA_320_NA_1_COD 0.0 0
#> 8 1963_1_5_1963_3_20_5_MAZ_44G9_1_NA_650_NA_1_COD 0.0 0
#> 9 1963_1_5_1963_3_20_5_MAZ_44G9_10_NA_390_NA_1_COD 0.0 0
#> 10 1963_1_5_1963_3_20_5_MAZ_44G9_11_NA_210_NA_1_COD 0.0 0
#> 11 1963_1_5_1963_3_20_5_MAZ_44G9_13_NA_220_NA_1_COD 0.0 0
#> 12 1963_1_5_1963_3_20_5_MAZ_44G9_14_NA_240_NA_1_COD 0.0 0
#> 13 1963_1_5_1963_3_20_5_MAZ_44G9_16_NA_230_NA_1_COD 0.0 0
#> 14 1963_1_5_1963_3_20_5_MAZ_44G9_17_NA_230_NA_1_COD 0.0 0
#> 15 1963_1_5_1963_3_20_5_MAZ_44G9_18_NA_160_NA_1_COD 0.0 0
#> 16 1963_1_5_1963_3_20_5_MAZ_44G9_20_NA_310_NA_1_COD 0.0 0
#> 17 1963_1_5_1963_3_20_5_MAZ_44G9_21_NA_240_NA_1_COD 0.0 0
#> 18 1963_1_5_1963_3_20_5_MAZ_44G9_22_NA_310_NA_1_COD 0.0 0
#> 19 1963_1_5_1963_3_20_5_MAZ_44G9_23_NA_410_NA_3_COD 0.0 0
#> 20 1963_1_5_1963_3_20_5_MAZ_44G9_27_NA_260_NA_1_COD 0.0 0
#> clupeidae_tot sprattus_sprattus_tot clupea_harengus_tot gadus_morhua_tot
#> 1 0 0 0.00 0
#> 2 0 0 0.00 0
#> 3 0 0 0.00 0
#> 4 0 0 0.00 0
#> 5 0 0 0.00 0
#> 6 0 0 0.00 0
#> 7 0 0 0.00 0
#> 8 0 0 48.77 0
#> 9 0 0 0.00 0
#> 10 0 0 0.00 0
#> 11 0 0 0.00 0
#> 12 0 0 0.00 0
#> 13 0 0 0.00 0
#> 14 0 0 0.00 0
#> 15 0 0 0.00 0
#> 16 0 0 0.00 0
#> 17 0 0 0.00 0
#> 18 0 0 0.00 0
#> 19 0 0 4.17 0
#> 20 0 0 0.00 0
#> gadiformes_tot gobiidae_tot mysidae_tot non_bio_tot other_crustacea_tot
#> 1 0 0 0 0 0.00
#> 2 0 0 0 0 0.00
#> 3 0 0 0 0 0.00
#> 4 0 0 0 0 0.00
#> 5 0 0 0 0 0.00
#> 6 0 0 0 0 0.00
#> 7 0 0 0 0 3.34
#> 8 0 0 0 0 0.00
#> 9 0 0 0 0 0.00
#> 10 0 0 0 0 0.10
#> 11 0 0 0 0 0.85
#> 12 0 0 0 0 0.97
#> 13 0 0 0 0 0.82
#> 14 0 0 0 0 0.32
#> 15 0 0 0 0 0.87
#> 16 0 0 0 0 1.00
#> 17 0 0 0 0 0.00
#> 18 0 0 0 0 0.00
#> 19 0 0 0 0 0.00
#> 20 0 0 0 0 0.00
#> other_tot other_pisces_tot platichthys_flesus_tot polychaeta_tot
#> 1 0 0 0 0.0
#> 2 0 0 0 1.2
#> 3 0 0 0 0.0
#> 4 0 0 0 0.0
#> 5 0 0 0 0.0
#> 6 0 0 0 0.0
#> 7 0 0 0 0.0
#> 8 0 0 0 0.0
#> 9 0 0 0 0.0
#> 10 0 0 0 0.0
#> 11 0 0 0 0.0
#> 12 0 0 0 0.0
#> 13 0 0 0 0.0
#> 14 0 0 0 0.0
#> 15 0 0 0 0.0
#> 16 0 0 0 0.0
#> 17 0 0 0 0.0
#> 18 0 0 0 0.0
#> 19 0 0 0 0.0
#> 20 0 0 0 0.0
#> saduria_entomon_tot year quarter cruise hn
#> 1 18.27 1963 1 No name 3
#> 2 5.62 1963 1 No name 3
#> 3 3.09 1963 1 No name 3
#> 4 12.60 1963 1 No name 3
#> 5 2.32 1963 1 No name 3
#> 6 1.50 1963 1 No name 3
#> 7 1.76 1963 1 No name 3
#> 8 0.00 1963 1 No name 5
#> 9 5.90 1963 1 No name 5
#> 10 0.67 1963 1 No name 5
#> 11 1.54 1963 1 No name 5
#> 12 0.00 1963 1 No name 5
#> 13 0.00 1963 1 No name 5
#> 14 0.00 1963 1 No name 5
#> 15 0.00 1963 1 No name 5
#> 16 0.45 1963 1 No name 5
#> 17 4.60 1963 1 No name 5
#> 18 3.37 1963 1 No name 5
#> 19 0.00 1963 1 No name 5
#> 20 2.85 1963 1 No name 5
#> sample predator_spec x y lat
#> 1 1963_3_18_3_MAZ_41G8_22_NA_500_NA_1 COD 716.8245 6239.414 56.25
#> 2 1963_3_18_3_MAZ_41G8_39_NA_450_NA_1 COD 716.8245 6239.414 56.25
#> 3 1963_3_18_3_MAZ_41G8_41_NA_560_NA_1 COD 716.8245 6239.414 56.25
#> 4 1963_3_18_3_MAZ_41G8_42_NA_390_NA_1 COD 716.8245 6239.414 56.25
#> 5 1963_3_18_3_MAZ_41G8_45_NA_270_NA_1 COD 716.8245 6239.414 56.25
#> 6 1963_3_18_3_MAZ_41G8_47_NA_270_NA_1 COD 716.8245 6239.414 56.25
#> 7 1963_3_18_3_MAZ_41G8_51_NA_320_NA_1 COD 716.8245 6239.414 56.25
#> 8 1963_3_20_5_MAZ_44G9_1_NA_650_NA_1 COD 767.7241 6409.776 57.75
#> 9 1963_3_20_5_MAZ_44G9_10_NA_390_NA_1 COD 767.7241 6409.776 57.75
#> 10 1963_3_20_5_MAZ_44G9_11_NA_210_NA_1 COD 767.7241 6409.776 57.75
#> 11 1963_3_20_5_MAZ_44G9_13_NA_220_NA_1 COD 767.7241 6409.776 57.75
#> 12 1963_3_20_5_MAZ_44G9_14_NA_240_NA_1 COD 767.7241 6409.776 57.75
#> 13 1963_3_20_5_MAZ_44G9_16_NA_230_NA_1 COD 767.7241 6409.776 57.75
#> 14 1963_3_20_5_MAZ_44G9_17_NA_230_NA_1 COD 767.7241 6409.776 57.75
#> 15 1963_3_20_5_MAZ_44G9_18_NA_160_NA_1 COD 767.7241 6409.776 57.75
#> 16 1963_3_20_5_MAZ_44G9_20_NA_310_NA_1 COD 767.7241 6409.776 57.75
#> 17 1963_3_20_5_MAZ_44G9_21_NA_240_NA_1 COD 767.7241 6409.776 57.75
#> 18 1963_3_20_5_MAZ_44G9_22_NA_310_NA_1 COD 767.7241 6409.776 57.75
#> 19 1963_3_20_5_MAZ_44G9_23_NA_410_NA_3 COD 767.7241 6409.776 57.75
#> 20 1963_3_20_5_MAZ_44G9_27_NA_260_NA_1 COD 767.7241 6409.776 57.75
#> long ices_rect pred_size_mm pred_weight_g source tot_prey_biom pred_cm
#> 1 18.5 41G8 500 NA Kevin's MSC 18.27 50
#> 2 18.5 41G8 450 NA Kevin's MSC 8.02 45
#> 3 18.5 41G8 560 NA Kevin's MSC 3.09 56
#> 4 18.5 41G8 390 NA Kevin's MSC 12.60 39
#> 5 18.5 41G8 270 NA Kevin's MSC 2.32 27
#> 6 18.5 41G8 270 NA Kevin's MSC 1.50 27
#> 7 18.5 41G8 320 NA Kevin's MSC 5.10 32
#> 8 19.5 44G9 650 NA Kevin's MSC 48.77 65
#> 9 19.5 44G9 390 NA Kevin's MSC 5.90 39
#> 10 19.5 44G9 210 NA Kevin's MSC 0.77 21
#> 11 19.5 44G9 220 NA Kevin's MSC 2.39 22
#> 12 19.5 44G9 240 NA Kevin's MSC 0.97 24
#> 13 19.5 44G9 230 NA Kevin's MSC 0.82 23
#> 14 19.5 44G9 230 NA Kevin's MSC 0.32 23
#> 15 19.5 44G9 160 NA Kevin's MSC 0.87 16
#> 16 19.5 44G9 310 NA Kevin's MSC 1.45 31
#> 17 19.5 44G9 240 NA Kevin's MSC 4.60 24
#> 18 19.5 44G9 310 NA Kevin's MSC 3.37 31
#> 19 19.5 44G9 410 NA Kevin's MSC 4.17 41
#> 20 19.5 44G9 260 NA Kevin's MSC 2.85 26
#> pred_cm_class
#> 1 (40,50]
#> 2 (40,50]
#> 3 (50,200]
#> 4 (30,40]
#> 5 (20,30]
#> 6 (20,30]
#> 7 (30,40]
#> 8 (50,200]
#> 9 (30,40]
#> 10 (20,30]
#> 11 (20,30]
#> 12 (20,30]
#> 13 (20,30]
#> 14 (20,30]
#> 15 (6,20]
#> 16 (30,40]
#> 17 (20,30]
#> 18 (30,40]
#> 19 (40,50]
#> 20 (20,30]
unique(cod$pred_cm_class)
#> [1] (40,50] (50,200] (30,40] (20,30] (6,20] <NA> (0,6]
#> Levels: (0,6] (6,20] (20,30] (30,40] (40,50] (50,200]
# First read the density models and predict the values
mcod2 <- readRDS("output/mcod2.rds")
mfle2 <- readRDS("output/mfle2.rds")
# Make predictions onto the stomach data defined earlier (1 row per stomach and prey?)
# Add depth and rename coordinates to match the data used for fitting the density model
# Read the tifs
west <- raster("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/cod_condition/data/depth_geo_tif/D5_2018_rgb-1.tif")
#plot(west)
east <- raster("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/cod_condition/data/depth_geo_tif/D6_2018_rgb-1.tif")
# plot(east)
dep_rast <- raster::merge(west, east)
colnames(cod)
#> [1] "unique_pred_id" "amphipoda_tot" "bivalvia_tot"
#> [4] "clupeidae_tot" "sprattus_sprattus_tot" "clupea_harengus_tot"
#> [7] "gadus_morhua_tot" "gadiformes_tot" "gobiidae_tot"
#> [10] "mysidae_tot" "non_bio_tot" "other_crustacea_tot"
#> [13] "other_tot" "other_pisces_tot" "platichthys_flesus_tot"
#> [16] "polychaeta_tot" "saduria_entomon_tot" "year"
#> [19] "quarter" "cruise" "hn"
#> [22] "sample" "predator_spec" "x"
#> [25] "y" "lat" "long"
#> [28] "ices_rect" "pred_size_mm" "pred_weight_g"
#> [31] "source" "tot_prey_biom" "pred_cm"
#> [34] "pred_cm_class" "FR_tot" "prop_saduria"
#> [37] "prop_common" "common_tot" "time_period"
#> [40] "quarter_fact"
colnames(fle)
#> [1] "unique_pred_id" "amphipoda_tot" "bivalvia_tot"
#> [4] "clupeidae_tot" "sprattus_sprattus_tot" "clupea_harengus_tot"
#> [7] "gadus_morhua_tot" "gadiformes_tot" "gobiidae_tot"
#> [10] "mysidae_tot" "non_bio_tot" "other_crustacea_tot"
#> [13] "other_tot" "other_pisces_tot" "platichthys_flesus_tot"
#> [16] "polychaeta_tot" "saduria_entomon_tot" "year"
#> [19] "quarter" "cruise" "hn"
#> [22] "sample" "predator_spec" "x"
#> [25] "y" "lat" "long"
#> [28] "ices_rect" "pred_size_mm" "pred_weight_g"
#> [31] "source" "tot_prey_biom" "pred_cm"
#> [34] "pred_cm_class" "FR_tot" "prop_saduria"
#> [37] "prop_common" "common_tot"
cod$depth <- extract(dep_rast, cod[, 27:26])
fle$depth <- extract(dep_rast, fle[, 27:26])
# Convert to depth (instead of elevation)
ggplot(cod, aes(depth)) + geom_histogram()
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
cod$depth <- (cod$depth - max(drop_na(cod)$depth)) *-1
#> drop_na: removed 1,658 rows (19%), 7,013 rows remaining
ggplot(cod, aes(depth)) + geom_histogram()
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
fle$depth <- (fle$depth - max(drop_na(fle)$depth)) *-1
#> drop_na: removed 556 rows (25%), 1,692 rows remaining
cod <- cod %>% rename("X" = "x", "Y" = "y")
#> rename: renamed 2 variables (X, Y)
cod <- cod %>% mutate(year = as.integer(year))
#> mutate: converted 'year' from double to integer (0 new NA)
fle <- fle %>% rename("X" = "x", "Y" = "y")
#> rename: renamed 2 variables (X, Y)
fle <- fle %>% mutate(year = as.integer(year))
#> mutate: converted 'year' from double to integer (0 new NA)
# Standardize depth to the DENSITY data, not this data set
density_dat <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/mdat_cpue.csv")
#>
#> ── Column specification ────────────────────────────────────────────────────────
#> cols(
#> density = col_double(),
#> year = col_double(),
#> lat = col_double(),
#> lon = col_double(),
#> quarter = col_double(),
#> haul.id = col_character(),
#> IDx = col_character(),
#> ices_rect = col_character(),
#> SD = col_double(),
#> density_fle = col_double(),
#> depth = col_double(),
#> oxy = col_double(),
#> temp = col_double(),
#> X = col_double(),
#> Y = col_double()
#> )
# Rename to match variables in the density model
cod <- cod %>% mutate(depth_sc = (depth - mean(density_dat$depth)) / sd(density_dat$depth))
#> mutate: new variable 'depth_sc' (double) with 106 unique values and 0% NA
fle <- fle %>% mutate(depth_sc = (depth - mean(density_dat$depth)) / sd(density_dat$depth))
#> mutate: new variable 'depth_sc' (double) with 81 unique values and 0% NA
cod2 <- cod %>% filter(year < 2020 & year > 1992) # Remove this after I've refitted the models with 2020!
#> filter: removed 2,029 rows (23%), 6,642 rows remaining
cod_old <- cod %>% filter(year < 1993) # Filter the old data so that we can rbind it later
#> filter: removed 6,914 rows (80%), 1,757 rows remaining
fle2 <- fle %>% filter(year < 2020) # Remove this after I've refitted the models with 2020!
#> filter: removed 481 rows (21%), 1,767 rows remaining
# Now predict using the density models
cod_stomach_predcod_density <- predict(mcod2, newdata = dplyr::select(cod2, depth_sc, year, X, Y))
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k`, specifing `k` now indicates
#> an upper limit on the basis functions, and the degree of wiggliness
#> is determined by the data.
cod_stomach_predfle_density <- predict(mfle2, newdata = dplyr::select(cod2, depth_sc, year, X, Y))
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k`, specifing `k` now indicates
#> an upper limit on the basis functions, and the degree of wiggliness
#> is determined by the data.
fle_stomach_predcod_density <- predict(mcod2, newdata = dplyr::select(fle2, depth_sc, year, X, Y))
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k`, specifing `k` now indicates
#> an upper limit on the basis functions, and the degree of wiggliness
#> is determined by the data.
fle_stomach_predfle_density <- predict(mfle2, newdata = dplyr::select(fle2, depth_sc, year, X, Y))
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k`, specifing `k` now indicates
#> an upper limit on the basis functions, and the degree of wiggliness
#> is determined by the data.
cod2$predcod_density <- exp(cod_stomach_predcod_density$est)
cod2$predfle_density <- exp(cod_stomach_predfle_density$est)
fle2$predcod_density <- exp(fle_stomach_predcod_density$est)
fle2$predfle_density <- exp(fle_stomach_predfle_density$est)
# Add scaled predicted densities to the stomach data
cod2 <- cod2 %>%
mutate(predcod_density_sc = (predcod_density - mean(predcod_density))/sd(predcod_density),
predfle_density_sc = (predfle_density - mean(predfle_density))/sd(predfle_density))
#> mutate: new variable 'predcod_density_sc' (double) with 441 unique values and 0% NA
#> new variable 'predfle_density_sc' (double) with 441 unique values and 0% NA
fle2 <- fle2 %>%
mutate(predcod_density_sc = (predcod_density - mean(predcod_density))/sd(predcod_density),
predfle_density_sc = (predfle_density - mean(predfle_density))/sd(predfle_density))
#> mutate: new variable 'predcod_density_sc' (double) with 192 unique values and 0% NA
#> new variable 'predfle_density_sc' (double) with 193 unique values and 0% NA
# Add the old data back using bind_rows so that density estimates get NA
cod2 <- bind_rows(cod2, cod_old)
colnames(cod)
#> [1] "unique_pred_id" "amphipoda_tot" "bivalvia_tot"
#> [4] "clupeidae_tot" "sprattus_sprattus_tot" "clupea_harengus_tot"
#> [7] "gadus_morhua_tot" "gadiformes_tot" "gobiidae_tot"
#> [10] "mysidae_tot" "non_bio_tot" "other_crustacea_tot"
#> [13] "other_tot" "other_pisces_tot" "platichthys_flesus_tot"
#> [16] "polychaeta_tot" "saduria_entomon_tot" "year"
#> [19] "quarter" "cruise" "hn"
#> [22] "sample" "predator_spec" "X"
#> [25] "Y" "lat" "long"
#> [28] "ices_rect" "pred_size_mm" "pred_weight_g"
#> [31] "source" "tot_prey_biom" "pred_cm"
#> [34] "pred_cm_class" "FR_tot" "prop_saduria"
#> [37] "prop_common" "common_tot" "time_period"
#> [40] "quarter_fact" "depth" "depth_sc"
cod_save <- cod2 %>% dplyr::select(FR_tot, saduria_entomon_tot, tot_prey_biom, common_tot, unique_pred_id,
year, quarter, time_period, quarter_fact, pred_weight_g, pred_cm,
predator_spec, predcod_density, predfle_density, predcod_density_sc,
predfle_density_sc, depth, X, Y, lat, long, ices_rect, cruise)
fle_save <- fle2 %>% dplyr::select(FR_tot, saduria_entomon_tot, tot_prey_biom, common_tot, unique_pred_id,
year, quarter, pred_weight_g, pred_cm,
predator_spec, predcod_density, predfle_density, predcod_density_sc,
predfle_density_sc, depth, X, Y, lat, long, ices_rect, cruise)
write.csv(cod_save, "data/cod_diet_analysis.csv")
write.csv(fle_save, "data/fle_diet_analysis.csv")
# How does stomach content scale with size? Probably a little!
ggplot(cod, aes(pred_cm, FR_tot, color = time_period)) +
geom_point(alpha = 0.4) +
facet_wrap(~time_period) +
stat_smooth(color = "black")
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
ggplot(fle, aes(pred_cm, FR_tot)) +
geom_point(alpha = 0.4) +
stat_smooth(color = "black")
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
# r plot cod stomach against density
# Common prey
p1 <- ggplot(cod2, aes(predcod_density_sc, prop_common)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth()
p2 <- ggplot(cod2, aes(predcod_density_sc, common_tot/pred_weight_g)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth()
p3 <- ggplot(cod2, aes(predfle_density_sc, prop_common)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth()
p4 <- ggplot(cod2, aes(predfle_density_sc, common_tot/pred_weight_g)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth()
p5 <- ggplot(cod2, aes(predfle_density_sc + predcod_density_sc, prop_common)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth()
p6 <- ggplot(cod2, aes(predfle_density_sc + predcod_density_sc, common_tot/pred_weight_g)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth()
(p1|p2) / (p3|p4) / (p5|p6) + plot_annotation(title = "Cod stomachs",
theme = theme(plot.title = element_text(size = 16)))
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> Warning: Removed 2467 rows containing non-finite values (stat_smooth).
#> Warning: Removed 2467 rows containing missing values (geom_point).
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> Warning: Removed 1757 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1757 rows containing missing values (geom_point).
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> Warning: Removed 2467 rows containing non-finite values (stat_smooth).
#> Warning: Removed 2467 rows containing missing values (geom_point).
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> Warning: Removed 1757 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1757 rows containing missing values (geom_point).
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> Warning: Removed 2467 rows containing non-finite values (stat_smooth).
#> Warning: Removed 2467 rows containing missing values (geom_point).
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> Warning: Removed 1757 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1757 rows containing missing values (geom_point).
# Saduria
p7 <- ggplot(cod2, aes(predcod_density_sc, prop_saduria)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth()
p8 <- ggplot(cod2, aes(predcod_density_sc, saduria_entomon_tot/pred_weight_g)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth()
p9 <- ggplot(cod2, aes(predfle_density_sc, prop_saduria)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth()
p10 <- ggplot(cod2, aes(predfle_density_sc, saduria_entomon_tot/pred_weight_g)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth()
p11 <- ggplot(cod2, aes(predfle_density_sc + predcod_density_sc, prop_saduria)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth()
p12 <- ggplot(cod2, aes(predfle_density_sc + predcod_density_sc, saduria_entomon_tot/pred_weight_g)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth()
(p7|p8) / (p9|p10) / (p11|p12) + plot_annotation(title = "Cod stomachs",
theme = theme(plot.title = element_text(size = 16)))
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> Warning: Removed 2467 rows containing non-finite values (stat_smooth).
#> Warning: Removed 2467 rows containing missing values (geom_point).
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> Warning: Removed 1757 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1757 rows containing missing values (geom_point).
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> Warning: Removed 2467 rows containing non-finite values (stat_smooth).
#> Warning: Removed 2467 rows containing missing values (geom_point).
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> Warning: Removed 1757 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1757 rows containing missing values (geom_point).
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> Warning: Removed 2467 rows containing non-finite values (stat_smooth).
#> Warning: Removed 2467 rows containing missing values (geom_point).
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> Warning: Removed 1757 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1757 rows containing missing values (geom_point).
# Total stomach content (FR)
p13 <- ggplot(cod2, aes(predcod_density_sc, FR_tot)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth()
p14 <- ggplot(cod2, aes(predfle_density_sc, FR_tot)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth()
p15 <- ggplot(cod2, aes(predfle_density_sc + predcod_density, FR_tot)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth()
(p13/p14/p15) + plot_annotation(title = "Cod stomachs", theme = theme(plot.title = element_text(size = 16)))
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> Warning: Removed 1757 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1757 rows containing missing values (geom_point).
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> Warning: Removed 1757 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1757 rows containing missing values (geom_point).
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> Warning: Removed 1757 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1757 rows containing missing values (geom_point).
# Closer inspection of some of the plots...
ggplot(cod2, aes(predfle_density_sc, prop_saduria)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(method = "lm") +
ggtitle("prop saduria vs fle density")
#> `geom_smooth()` using formula 'y ~ x'
#> Warning: Removed 2467 rows containing non-finite values (stat_smooth).
#> Warning: Removed 2467 rows containing missing values (geom_point).
ggplot(cod2, aes(predfle_density, tot_prey_biom/pred_weight_g)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth() +
ggtitle("g/g tot. stomach vs fle density")
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> Warning: Removed 1757 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1757 rows containing missing values (geom_point).
ggplot(cod2, aes(predfle_density + predcod_density, tot_prey_biom/pred_weight_g)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth() +
ggtitle("g/g tot. stomach vs fle + coddensity")
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> Warning: Removed 1757 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1757 rows containing missing values (geom_point).
## Plot flounder stomach against density
# Common prey
p1 <- ggplot(fle2, aes(predcod_density_sc, prop_common)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth()
p2 <- ggplot(fle2, aes(predcod_density_sc, common_tot/pred_weight_g)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth()
p3 <- ggplot(fle2, aes(predfle_density_sc, prop_common)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth()
p4 <- ggplot(fle2, aes(predfle_density_sc, common_tot/pred_weight_g)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth()
p5 <- ggplot(fle2, aes(predfle_density_sc + predcod_density_sc, prop_common)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth()
p6 <- ggplot(fle2, aes(predfle_density_sc + predcod_density_sc, common_tot/pred_weight_g)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth()
(p1|p2) / (p3|p4) / (p5|p6) + plot_annotation(title = "Flounder stomachs",
theme = theme(plot.title = element_text(size = 16)))
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> Warning: Removed 516 rows containing non-finite values (stat_smooth).
#> Warning: Removed 516 rows containing missing values (geom_point).
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> Warning: Removed 516 rows containing non-finite values (stat_smooth).
#> Warning: Removed 516 rows containing missing values (geom_point).
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> Warning: Removed 516 rows containing non-finite values (stat_smooth).
#> Warning: Removed 516 rows containing missing values (geom_point).
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
# Saduria
p7 <- ggplot(fle2, aes(predcod_density_sc, prop_saduria)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth()
p8 <- ggplot(fle2, aes(predcod_density_sc, saduria_entomon_tot/pred_weight_g)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth()
p9 <- ggplot(fle2, aes(predfle_density_sc, prop_saduria)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth()
p10 <- ggplot(fle2, aes(predfle_density_sc, saduria_entomon_tot/pred_weight_g)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth()
p11 <- ggplot(fle2, aes(predfle_density_sc + predcod_density_sc, prop_saduria)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth()
p12 <- ggplot(fle2, aes(predfle_density_sc + predcod_density_sc, saduria_entomon_tot/pred_weight_g)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth()
(p7|p8) / (p9|p10) / (p11|p12) + plot_annotation(title = "Flounder stomachs",
theme = theme(plot.title = element_text(size = 16)))
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> Warning: Removed 516 rows containing non-finite values (stat_smooth).
#> Warning: Removed 516 rows containing missing values (geom_point).
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> Warning: Removed 516 rows containing non-finite values (stat_smooth).
#> Warning: Removed 516 rows containing missing values (geom_point).
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> Warning: Removed 516 rows containing non-finite values (stat_smooth).
#> Warning: Removed 516 rows containing missing values (geom_point).
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
# FR
p13 <- ggplot(fle2, aes(predcod_density_sc, FR_tot)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth()
p14 <- ggplot(fle2, aes(predfle_density_sc, FR_tot)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth()
p15 <- ggplot(fle2, aes(predfle_density_sc + predcod_density, FR_tot)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth()
(p13/p14/p15) + plot_annotation(title = "Flounder stomachs", theme = theme(plot.title = element_text(size = 16)))
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
# Does it seem like proportions of the diet changes with density, but not the total amount of food?
# Check saduria and flounder for instance.
pp7 <- ggplot(fle2, aes(predcod_density_sc, prop_saduria)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(method = "lm")
pp8 <- ggplot(filter(fle2, (saduria_entomon_tot/pred_weight_g) < 0.1),
aes(predcod_density_sc, saduria_entomon_tot/pred_weight_g)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(method = "lm")
#> filter: no rows removed
pp13 <- ggplot(filter(fle2, (saduria_entomon_tot/pred_weight_g) < 0.1),
aes(predcod_density_sc, FR_tot)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(method = "lm")
#> filter: no rows removed
ppfle <- pp7/pp8/pp13 +
plot_layout(ncol = 3) +
plot_annotation(title = "Flounder stomachs", theme = theme(plot.title = element_text(size = 16)))
# And that cod density is unaffected (or less at least)
p7 <- ggplot(filter(cod2, (saduria_entomon_tot/pred_weight_g) < 0.02 & predfle_density_sc < 2),
aes(predfle_density_sc, prop_saduria)) + # Crop it up a little!
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(method = "lm")
#> filter: removed 2,036 rows (24%), 6,363 rows remaining
p8 <- ggplot(filter(cod2, (saduria_entomon_tot/pred_weight_g) < 0.02 & predfle_density_sc < 2),
aes(predfle_density_sc, saduria_entomon_tot/pred_weight_g)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(method = "lm")
#> filter: removed 2,036 rows (24%), 6,363 rows remaining
p13 <- ggplot(filter(cod2, (saduria_entomon_tot/pred_weight_g) < 0.02 & predfle_density_sc < 2 & FR_tot < 0.5),
aes(predfle_density_sc, FR_tot)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(method = "lm")
#> filter: removed 2,036 rows (24%), 6,363 rows remaining
ppcod <- p7/p8/p13 +
plot_layout(ncol = 3) +
plot_annotation(title = "Cod stomachs", theme = theme(plot.title = element_text(size = 16)))
ppfle / ppcod
#> `geom_smooth()` using formula 'y ~ x'
#> Warning: Removed 516 rows containing non-finite values (stat_smooth).
#> Warning: Removed 516 rows containing missing values (geom_point).
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> Warning: Removed 690 rows containing non-finite values (stat_smooth).
#> Warning: Removed 690 rows containing missing values (geom_point).
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
# What does it mean? Yes, the proportion *may* be negatively affected by densities,
# but the food content is not, which I guess would matter for condition
# Now, group by survey and Quarter, select only the main groups (as Haase) and see if it resembles the
# previous results.
## Reproduce Haase figures with biomass proportions
# Try to reproduce Haase et al!
# Aggregate across predator individuals and convert to long data for plotting
cod2 <- cod %>%
drop_na(pred_cm_class) %>% # This is because not all predators have length
filter(cruise %in% c("BITS", "BITS-1", "BITS-2")) %>%
filter(year %in% c(2015, 2016, 2017)) %>%
filter(source == "Kevin's MSC") %>%
filter(!pred_cm_class == "(0,6]") %>%
group_by(quarter, pred_cm_class) %>%
summarise(sprattus_sprattus_tot_all = sum(sprattus_sprattus_tot),
saduria_entomon_tot_all = sum(saduria_entomon_tot),
mysidae_tot_all = sum(mysidae_tot),
clupea_harengus_tot_all = sum(clupea_harengus_tot),
sprattus_sprattus_tot_all = sum(sprattus_sprattus_tot),
gobiidae_tot_all = sum(gobiidae_tot),
gadiformes_tot_all = sum(gadiformes_tot),
other_pisces_tot_all = sum(other_pisces_tot),
polychaeta_tot_all = sum(polychaeta_tot)) %>%
ungroup() %>%
pivot_longer(3:10) %>%
group_by(quarter, pred_cm_class) %>%
mutate(percent = value/sum(value)) %>%
ungroup()
#> drop_na: no rows removed
#> filter: removed 3,358 rows (39%), 5,313 rows remaining
#> filter: removed 3,191 rows (60%), 2,122 rows remaining
#> filter: removed 283 rows (13%), 1,839 rows remaining
#> filter: removed 16 rows (1%), 1,823 rows remaining
#> group_by: 2 grouping variables (quarter, pred_cm_class)
#> summarise: now 10 rows and 10 columns, one group variable remaining (quarter)
#> ungroup: no grouping variables
#> pivot_longer: reorganized (sprattus_sprattus_tot_all, saduria_entomon_tot_all, mysidae_tot_all, clupea_harengus_tot_all, gobiidae_tot_all, …) into (name, value) [was 10x10, now 80x4]
#> group_by: 2 grouping variables (quarter, pred_cm_class)
#> mutate (grouped): new variable 'percent' (double) with 64 unique values and 0% NA
#> ungroup: no grouping variables
# So much herring for small cod?
# A 20 cm cod weighs approx 100 g. Can it have 20g of herring in the stomach?
#small_cod <-
d %>%
filter(Pred.size.mm < 200 & Prey.sp.code %in% c("clupea harengus", "Clupea harengus")) %>%
dplyr::select(Prey.weight, Prey.sp.code, Unique.pred.id, Pred.size.mm)
#> filter: removed 25,589 rows (>99%), 46 rows remaining
#> # A tibble: 46 x 4
#> Prey.weight Prey.sp.code Unique.pred.id Pred.size.mm
#> <dbl> <chr> <chr> <dbl>
#> 1 20.0 Clupea hareng… 2016_1_73_735_COD 180
#> 2 18.4 Clupea hareng… 1977_1_16_1977_1_30_16_UNKN_41G8_39_… 0
#> 3 8.02 Clupea hareng… 1977_1_16_1977_1_30_16_UNKN_41G8_67_… 0
#> 4 0.44 Clupea hareng… 1977_1_16_1977_1_30_16_UNKN_41G8_68_… 0
#> 5 7.25 Clupea hareng… 1977_1_9_1977_3_15_9_UNKN_41G8_25_NA… 0
#> 6 26.2 Clupea hareng… 1977_1_9_1977_3_15_9_UNKN_41G8_26_NA… 0
#> 7 6.17 Clupea hareng… 1977_1_9_1977_3_15_9_UNKN_41G8_30_NA… 0
#> 8 18.9 Clupea hareng… 1977_1_9_1977_3_15_9_UNKN_41G8_31_NA… 0
#> 9 16.4 Clupea hareng… 1977_1_9_1977_3_15_9_UNKN_41G8_38_NA… 0
#> 10 44.2 Clupea hareng… 1977_1_9_1977_3_15_9_UNKN_41G8_41_NA… 0
#> # … with 36 more rows
d_wide4 %>% filter(unique_pred_id == "2016_1_73_735_COD") %>% as.data.frame()
#> filter: removed 11,188 rows (>99%), one row remaining
#> unique_pred_id amphipoda_tot bivalvia_tot clupeidae_tot
#> 1 2016_1_73_735_COD 0.05 0 0
#> sprattus_sprattus_tot clupea_harengus_tot gadus_morhua_tot gadiformes_tot
#> 1 0 19.99 0 0
#> gobiidae_tot mysidae_tot non_bio_tot other_crustacea_tot other_tot
#> 1 0 0 0 0 0
#> other_pisces_tot platichthys_flesus_tot polychaeta_tot saduria_entomon_tot
#> 1 0 0 0.05 0
#> year quarter cruise hn sample predator_spec x y lat long
#> 1 2016 1 BITS 73 735 COD 613.9639 6271.507 56.574 16.855
#> ices_rect pred_size_mm pred_weight_g source tot_prey_biom
#> 1 42G6 180 46 Kevin's MSC 20.09
d %>% filter(Unique.pred.id == "2016_1_73_735_COD") %>% as.data.frame()
#> filter: removed 25,633 rows (>99%), 2 rows remaining
#> SubDiv Year Month Day N.with.food N.regurgitated N.skeletal N.empty HN Sample
#> 1 27 2016 2 28 1 NA NA NA 73 735
#> 2 27 2016 2 28 NA NA 1 NA 73 735
#> Length.code Prey.sp.code Prey.size Prey.weight Prey.nr Stage.digestion
#> 1 Lt Clupea harengus 150 19.99 1 1
#> 2 <NA> Bylgides sarsi NA 0.05 NA 2
#> Comment Parasites.in.stomach Quarter Coun. Ship Cruise Pred.size.mm
#> 1 <NA> <NA> 1 SWE NA BITS 180
#> 2 <NA> <NA> 1 SWE NA BITS 180
#> Pred.weight.g Predator.gutted.weight Predator.code Sex Maturity Gall.content
#> 1 46 42 COD M 2 4
#> 2 46 42 COD M 2 4
#> Q.year Age Date Index Lat Long Depth.catch Ices.rect source
#> 1 1_2016 2 <NA> OXBH.2016.73 56.574 16.855 63 42G6 Kevin's MSC
#> 2 1_2016 2 <NA> OXBH.2016.73 56.574 16.855 63 42G6 Kevin's MSC
#> Number Sample.type transect Depthstep Gonad.weight.roundfish
#> 1 NA NA NA NA NA
#> 2 NA NA NA NA NA
#> Liver.weight.roundfish Stomach.content Perc.stomac.content comments Processed
#> 1 NA NA NA NA NA
#> 2 NA NA NA NA NA
#> Unique.pred.id coord_from_ices_rect_centre X Y
#> 1 2016_1_73_735_COD N 613.9639 6271.507
#> 2 2016_1_73_735_COD N 613.9639 6271.507
# Check sample sizes per size class (numbers below bars in Kevins figures)
cod %>%
drop_na(pred_cm_class) %>% # This is because not all predators have length
filter(cruise %in% c("BITS", "BITS-1", "BITS-2")) %>%
filter(year %in% c(2015, 2016, 2017)) %>%
filter(source == "Kevin's MSC") %>%
filter(!pred_cm_class == "(0,6]") %>%
group_by(quarter, pred_cm_class) %>%
summarise(n = n())
#> drop_na: no rows removed
#> filter: removed 3,358 rows (39%), 5,313 rows remaining
#> filter: removed 3,191 rows (60%), 2,122 rows remaining
#> filter: removed 283 rows (13%), 1,839 rows remaining
#> filter: removed 16 rows (1%), 1,823 rows remaining
#> group_by: 2 grouping variables (quarter, pred_cm_class)
#> summarise: now 10 rows and 3 columns, one group variable remaining (quarter)
#> # A tibble: 10 x 3
#> # Groups: quarter [2]
#> quarter pred_cm_class n
#> <dbl> <fct> <int>
#> 1 1 (6,20] 136
#> 2 1 (20,30] 275
#> 3 1 (30,40] 303
#> 4 1 (40,50] 194
#> 5 1 (50,200] 120
#> 6 4 (6,20] 128
#> 7 4 (20,30] 266
#> 8 4 (30,40] 236
#> 9 4 (40,50] 132
#> 10 4 (50,200] 33
# I seem to have more data...
# Check Kevin's raw data...
kev <- read.csv("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/bentfish/data/stomach_data/Master_Kevin_updated.csv", sep = ";")
head(kev)
#> X coun. ship cruise sample.typ sd ices rect transect date qyear
#> 1 1 SWE DANS BITS trawl SD26-28 43G8 4363 23.11.2015 4_2015
#> 2 2 SWE DANS BITS trawl SD26-28 43G8 4363 23.11.2015 4_2015
#> 3 3 SWE DANS BITS trawl SD26-28 43G8 4363 23.11.2015 4_2015
#> 4 4 SWE DANS BITS trawl SD26-28 43G8 4363 23.11.2015 4_2015
#> 5 5 SWE DANS BITS trawl SD26-28 43G8 4363 23.11.2015 4_2015
#> 6 6 SWE DANS BITS trawl SD26-28 43G8 4363 23.11.2015 4_2015
#> year month day quarter station.haul index lat long depth
#> 1 2015 11 23 4 19 OXBH.2015.19 57,005 18,806 72
#> 2 2015 11 23 4 19 OXBH.2015.19 57,005 18,806 72
#> 3 2015 11 23 4 19 OXBH.2015.19 57,005 18,806 72
#> 4 2015 11 23 4 19 OXBH.2015.19 57,005 18,806 72
#> 5 2015 11 23 4 19 OXBH.2015.19 57,005 18,806 72
#> 6 2015 11 23 4 19 OXBH.2015.19 57,005 18,806 72
#> depthstep sample pred.code pred.size length lengthclass pred.weight
#> 1 >30m 406 COD 190 19 Lc6-20 70
#> 2 >30m 385 COD 320 32 Lc31-40 246
#> 3 >30m 389 COD 330 33 Lc31-40 330
#> 4 >30m 396 COD 300 30 Lc21-30 193
#> 5 >30m 381 COD 310 31 Lc31-40 244
#> 6 >30m 400 COD 250 25 Lc21-30 136
#> pred.gutted.weight age sex maturity gonad.weight liver.weight gall.bladder
#> 1 60 1 M 2 NA 5 1
#> 2 213 3 M 8 NA 9 1
#> 3 499 3 M 2 NA 27 4
#> 4 169 2 M 3 NA 10 4
#> 5 220 2 M 3 NA 11 4
#> 6 115 1 M 3 NA 8 2
#> specimen.note n.full n.regurgitated n.skeletal n.empty stomachcontent
#> 1 1 NA NA NA 0,87
#> 2 1 NA NA NA 1,78
#> 3 NA 1 NA NA 0
#> 4 1 NA NA NA 0,42
#> 5 1 NA NA NA 0,06
#> 6 1 NA NA NA 0,8
#> perc.stomachcontent length.code prey.sp.code prey.size prey.weight prey.nr
#> 1 1,24 Mysis mixta NA 0,87 19
#> 2 0,72 Ls Clupea harengus 85 1,78 1
#> 3 0 NA 0
#> 4 0,22 Mysis mixta NA 0,42 10
#> 5 0,02 Mysis mixta NA 0,06 2
#> 6 0,59 Mysis mixta NA 0,8 19
#> stage.digestion comment parasites.in.stomach processed
#> 1 1 original
#> 2 1 without head original
#> 3 NA fresh scales sure
#> 4 1 original
#> 5 1 original
#> 6 1 original
kev$prey.weight <- as.numeric(gsub(",", ".", gsub("\\.", "", kev$prey.weight)))
kev %>%
filter(lengthclass == "Lc6-20" & pred.code == "COD" & year %in% c(2015, 2016, 2017) & quarter == 1) %>%
drop_na(prey.sp.code) %>%
drop_na(prey.weight) %>%
group_by(prey.sp.code) %>%
summarise(tot_weight_by_prey = sum(prey.weight)) %>%
arrange(desc(tot_weight_by_prey)) %>%
as.data.frame()
#> filter: removed 15,869 rows (99%), 223 rows remaining
#> drop_na: no rows removed
#> drop_na: no rows removed
#> group_by: one grouping variable (prey.sp.code)
#> summarise: now 23 rows and 2 columns, ungrouped
#> prey.sp.code tot_weight_by_prey
#> 1 Clupea harengus 19.99
#> 2 Mysis mixta 6.52
#> 3 Bylgides sarsi 4.52
#> 4 Saduria entomon 3.75
#> 5 Gobiidae 3.53
#> 6 Clupeidae 1.29
#> 7 Diastylis rathkei 1.11
#> 8 Crangon crangon 0.58
#> 9 Priapulida 0.40
#> 10 Gammarus sp. 0.28
#> 11 Mysidae 0.22
#> 12 Halicryptus spinulosus 0.17
#> 13 Cumacea 0.15
#> 14 Limecola balthica 0.11
#> 15 Bivalvia 0.10
#> 16 Pontoporeia femorata 0.07
#> 17 Remains 0.06
#> 18 Amphipoda 0.01
#> 19 Monoporeia affinis 0.01
#> 20 Mytilus sp. 0.01
#> 21 Neomysis integer 0.01
#> 22 0.00
#> 23 Copepoda 0.00
# OK so probably Kevin made some additional filter, because I also have a larger sample size